/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
    \file integrator.hpp
    \brief Gaussian integration class.

    This class implements the Gaussian integration on the GLL nodes on a given
    rectangular dominium. An affine transformation is used to transoform the
    domain in the square \f$ [-1,1] \times [-1,1] \f$.
*/

# ifndef INTEGRATOR_HPP
# define INTEGRATOR_HPP 1

#include<cstdlib>
#include<cassert>

#include "point.hpp"
#include "poly.hpp"
#include "functions.hpp"
#include "typedefs.hpp"

namespace spectral {

/**
    \class Integrator
    \brief GLL Integration Class

    This class implements the Gaussian integration on the GLL nodes on a given
    rectangular dominium.
*/

/**
    \fn Integrator::Integrator(int nnodes, double delta, double eps, int nmax, int order)
    \brief Integrator constructor

    This will initializate the Integrator with \e nnodes nodes
    over a rectangular dominium in the form:
    \f[ [A_x,B_x] \times [A_y,B_y]\f]
    delta is used as a stop paramether on the increment for the newton method
    eps is used as a residual stop paramether.

    The lagrange basis is prebuild and then stored.

    \e nmax is the maximum number of iteration for the newton method, while
    computing the roots of the Legendre polynomials.

    \e order is the order of the problem in the week form.
    For example it's 1 for the standard advection-diffusion problem.

    \e nnodes must be at least 4 or our method for the generation
    of the lagrangian basis could fail.
*/

/**
    \fn Integrator::Integrator(Integrator const & integr)
    \brief Integrator copy constructor.

*/

/**
    \fn Integrator::~Integrator()
    \brief Integrator destructor
*/

/**
    \fn double Integrator::operator()(T const & f, Point<double> const & A,
                                      Point<double> const & B) const
    \brief Integrates

    This will integrate the function \e f over the given domain, using
    the gaussian integration over the GLL nodes.
    The template parameter \e T must have an overloaded operator "()" with
    signature:
    \code
    double operator()(double,double)
    \endcode
*/

/**
    \fn double Integrator::integrate(double const * v,
                                     Point<double> const & A,
                                     Point<double> const & B) const
    \brief Integrates

    This will integrate the vector \e v over the given domain, using
    the gaussian integration over the GLL nodes.
    \e v must be a vector containing the values of the function on the domain.

    This method is to be used to compute the residual L2 norm of a vector.
    The functor approach is to be prefered for other computations.
*/

/**
    \fn double Integrator::operator()(T const & f,
				      double const & a,
				      double const & b,
				      double const & xi,
				      unsigned int const & border) const
    \brief Integrates

    This will integrate the 2d function \f$ f(x,y) \f$ over the given 1d
    dominium \f$ [a,b] \f$ on the border \c border, using the gaussian
    integration over the GLL nodes. The parameter \c xi must contain the
    value of the variable which is not integrated on the border \c border.
    In order to make this method work, the template parameter \c T must
    have an overloaded operator "()" with signature:
    \code
    double operator()(double,double)
    \endcode
*/

/**
    \fn TensorProd Integrator::lagrange(in i) const
    \brief Return \f$ l_i \f$

    This returns an element of the 1d lagrangian basis on the reference interval
    \f$ [-1,1] \f$
*/

/**
    \fn TensorProd Integrator::lagrange(int i, int j) const
    \brief Returns \f$ l_{ij} \f$

    This returns an element of the 2d lagrangian basis on the reference square
    \f$ [-1,1] \times [-1,1] \f$
*/

/**
    \fn double * Integrator::getnodes() const
    \brief Returns the integration nodes

    This returns a pointer to an array of double containig the values of the
    quadrature nodes.
*/

/**
    \fn TensorProd Integrator::Dx(TensorProd const & tp,double lenghtX) const
    \brief Returns the partial derivative of tp

    This returns the partial derivative with respect to x on the physical domain
    of tp. \e lengthX must be the length of the physical element and it is used
    to compute the correct weight for the derivative:
    \f[
	\frac{\partial f}{\partial x} = \frac{\partial f}{\partial \hat{x}}
					\frac{\partial \hat{x}}{\partial x}
    \f]
*/

/**
    \fn TensorProd Integrator::Dy(TensorProd const & tp,double lenghtY) const
    \brief Returns the partial derivative of tp

    This returns the partial derivative with respect to x of the physical domain
    of tp. \e lengthY must be the length of the physical element and it's used
    to compute the correct weight for the derivative:
    \f[
	\frac{\partial f}{\partial y} = \frac{\partial f}{\partial \hat{y}}
					\frac{\partial \hat{y}}{\partial y}
    \f]
*/

/**
   \fn double Integrator::getnode(int n) const
   \brief Returns the integrator node \e n

   This returns the value of the quadrature node \e n
*/

/**
    \fn int Integrator::getnumnodes() const
    \brief Number of nodes for the quadrature formula.

*/

class Integrator {
    public:
        Integrator(int nnodes,
                   double delta,
                   double eps,
                   int nmax=1000,
                   int order=1);

        Integrator(Integrator const & integr);

        ~Integrator();

        template<class T>
        double operator()(T const & f,
                          Point<double> const & A,
                          Point<double> const & B) const;

        double integrate(double const * v,
                         Point<double> const & A,
                         Point<double> const & B) const;

        template<class T>
        double operator()(T const & f,
			  double const & a,
			  double const & b,
			  double const & xi,
			  unsigned int const & border) const;

	TensorProd lagrange(int i) const;
        TensorProd lagrange(int i, int j) const;

        TensorProd Dx(TensorProd const & tp, double lengthX) const;
        TensorProd Dy(TensorProd const & tp, double lengthY) const;

        int getnumnodes() const;
        double getnode(int n) const;
        double * getnodes() const;
    private:
        double * _weightsX;  /**< \brief Cached weights for the x-axis */
        double * _weightsY;  /**< \brief Cached weights for the y-axis */

        int _nnodes;         /**< \brief Number of nodes */
        double * _nodes;     /**< \brief GLL Nodes on [-1,1] */
        double * _weights;   /**< \brief GLL Weights on [-1,1] */
        double _delta;       /**< \brief Incremental stop paramether for the newton method */
        double _eps;         /**< \brief Residual stop paramether for the newton method */
        int _nmax;           /**< \brief Maximum number of iteration for the newton method */
        int _order;          /**< \brief Maximum order of derivation */
        Poly ** _lagrange;   /**< \brief Lagrange basis */

        static int _compare(const void * a, const void * b); /**< \brief Used for the C sort */
};

// ====================================================
// End of proper header
// ====================================================

template<class T>
double Integrator::operator()(T const & f,
                              Point<double> const & A,
                              Point<double> const & B) const{
    double lengthX=B.get(0)-A.get(0);
    assert(lengthX>0);
    double lengthY=B.get(1)-A.get(1);
    assert(lengthY>0);

    // Precalculating weights and nodes is more efficient
    // it's O(n) instead of O(n^2):
    for(int i=0;i<_nnodes;++i) {
       _weightsX[i]=0.5*lengthX*_weights[i];
       _weightsY[i]=0.5*lengthY*_weights[i];
    }

    double sum(0.0);

    for(int i=0;i<_nnodes;++i) {
        for(int j=0;j<_nnodes;++j) {
            sum+=_weightsX[i]*_weightsY[j]*f(_nodes[i],_nodes[j]);
        }
    }

    return sum;
}

template<class T>
double Integrator::operator()(T const & f,
			      double const & a,
			      double const & b,
			      double const & xi,
			      unsigned int const & border) const {
    double length=std::abs(b-a);

    double sum(0.0);

    if (border%2) {
	for(int i=0;i<_nnodes;++i)
	    sum+=0.5*length*_weights[i]*f(xi,_nodes[i]);
    } else {
	for(int i=0;i<_nnodes;++i)
	    sum+=0.5*length*_weights[i]*f(_nodes[i],xi);
    }

    return sum;
}

} // namespace

#endif
